t2 = linspace(0, 360, 360); 

S0_values = zeros(1, length(t2));

S_in = [1; 1; 0; 0];

g = 0;
t1 = 0;
f = 90;

for i = 1:length(t2)
    rotated_polarizer = 0.5 * [1, cosd(2*g)*cosd(2*t1), cosd(2*g)*sind(2*t1), 0;
                               cosd(2*g)*cosd(2*t1), cosd(2*t1)^2 + sind(2*g)*sind(2*t1)^2, (1 - sind(2*g))*sind(2*t1)*cosd(2*t1), 0;
                               cosd(2*g)*sind(2*t1), (1 - sind(2*g))*sind(2*t1)*cosd(2*t1), sind(2*t1)^2 + sind(2*g)*cosd(2*t1)^2, 0;
                               0, 0, 0, sind(2*g)];
    
    
    rotated_retarder = [1, 0, 0, 0;
                        0, (cosd(2*t2(i)))^2 + cosd(f)*sind(2*t2(i))^2, (1 - cosd(f))*sind(2*t2(i))*cosd(2*t2(i)), -sind(f)*sind(2*t2(i));
                        0, (1 - cosd(f))*sind(2*t2(i))*cosd(2*t2(i)), sind(2*t2(i))^2 + cosd(f)*cosd(2*t2(i))^2, sind(f)*cosd(2*t2(i));
                        0, sind(f)*sind(2*t2(i)), -sind(f)*cosd(2*t2(i)), cosd(f)];

    
    S_out = rotated_polarizer * rotated_retarder * S_in;
    
    
    S0_values(i) = S_out(1);
end

figure;
plot(t2, S0_values);
title('S0 as a function of t2 (1/4retarder)');
xlabel('t2 (degrees)');
ylabel('S0 (intensity)');
grid on;
